function x_new = rk4_capacitive(x0, dt, b, r, c, n, lambda_C, I_Chaos, u_PC)
    % rk4_capacitive - 使用四阶 Runge-Kutta 方法对电容耦合模型进行单步更新
    %
    % 输入:
    %   y0       - 当前状态向量 [V1; iL1; V2; iL2] (列向量)
    %   dt       - 单步步长
    %   b, r, c  - 无量纲参数 (例如 b = E_j/V0, r = R_j/a, c = a^2 * C_j / L_j)
    %   n        - 与压电传感器相关的参数
    %   lambda_C - 电容耦合强度（例如论文中给定 lambda_C = 0.2）
    %   I_Chaos  - 混沌电流（无量纲）
    %   u_PC     - 外部声信号（压电电压）
    %
    % 输出:
    %   y_new    - 单步更新后的状态向量 [V1; iL1; V2; iL2]

    % 确保 y0 为列向量
    x0 = x0(:);
    
    % 计算四个 RK4 阶段
    k1 = dt * capacitive(x0, b, r, c, n, lambda_C, I_Chaos, u_PC);
    k2 = dt * capacitive(x0 + 0.5 * k1, b, r, c, n, lambda_C, I_Chaos, u_PC);
    k3 = dt * capacitive(x0 + 0.5 * k2, b, r, c, n, lambda_C, I_Chaos, u_PC);
    k4 = dt * capacitive(x0 + k3, b, r, c, n, lambda_C, I_Chaos, u_PC);
    
    % 更新状态：使用 RK4 公式
    x_new = x0 + (k1 + 2*k2 + 2*k3 + k4) / 6;
end

function dXdt = capacitive(X, b, r, c, n, lambda_C, I_Chaos, u_PC)
    % capacitive_model - 计算电容耦合模型的状态导数
    %
    % 输入:
    %   X        - 当前状态向量 [V1; iL1; V2; iL2]
    %   b, r, c  - 无量纲参数
    %   n        - 与压电传感器相关的参数
    %   lambda_C - 电容耦合强度
    %   I_Chaos  - 混沌电流（无量纲）
    %   u_PC     - 外部声信号（压电电压）
    %
    % 输出:
    %   dXdt     - 状态导数 [dx1; dy1; dx2; dy2]
    
    % 解构状态变量
    x1 = X(1);
    y1 = X(2);
    x2 = X(3);
    y2 = X(4);
    
    % 计算两个中间项：
    % part1: x1 部分的非线性项，包含 I_Chaos
    part1 = x1 - (1/3)*x1^3 - y1 + I_Chaos;
    % part2: x2 部分的非线性项，包含外部信号（用 n * u_PC 表示）
    part2 = x2*(1 - n) - (1/3)*x2^3 - y2 + n*u_PC;
    
    % 根据论文方程（8）计算电容耦合调节项
    dx1 = part1 + lambda_C * (part2 - part1);
    dy1 = c * (x1 + b - r*y1);
    dx2 = part2 - lambda_C * (part2 - part1);
    dy2 = c * (x2 + b - r*y2);
    
    dXdt = [dx1; dy1; dx2; dy2];
end
